Systems of Nonlinear Equations

Problem. Let us find a solution $ξ = (ξ_{1}, ξ_{2}, \dots, ξ_{n})$ of the system of $n$ equations

$\begin{aligned} f_{1} (x) & = 0, \\ f_{2} (x) & = 0, \\ ⋮ \\ f_{n} (x) & = 0, \end{aligned}$

and $n$ unknowns $x = (x_{1}, x_{2}, \dots, x_{n})$ . Denoting $f = (f_{1}, f_{2}, \dots, f_{n})^{T}$ , we can write this system as

$f (x) = 0.$

We shall describe the Newton's method and three quasi-Newton methods:

Broyden's method,
Davidon-Fletcher-Powell method, and
Broyden-Fletcher-Goldfarb-Schano method.

Given starting approximation $x^{(0)}$ , all methods generate sequence of points $x^{(n)}$ which, under certain condditions, converges towards the solution $ξ$ .

Remark. Description of the methods and examples are taken from the textbook Numerička matematika, poglavlje 4.4.

3.3 ms

Newton's method

Jacobian or Jacobi matrix of functions $f$ in the point $x$ is the matrix of first partial derivatives

$J (f, x) = (\begin{matrix} \frac{\partial f_{1} (x)}{\partial x_{1}} & \frac{\partial f_{1} (x)}{\partial x_{2}} & \dots & \frac{\partial f_{1} (x)}{\partial x_{n}} \\ \frac{\partial f_{2} (x)}{\partial x_{1}} & \frac{\partial f_{2} (x)}{\partial x_{2}} & \dots & \frac{\partial f_{2} (x)}{\partial x_{n}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \frac{\partial f_{n} (x)}{\partial x_{1}} & \frac{\partial f_{n} (x)}{\partial x_{2}} & \dots & \frac{\partial f_{n} (x)}{\partial x_{n}} \end{matrix}) .$

Given starting approximation $x^{(0)}$ , we compute the sequence of points

$x^{(k + 1)} = x^{(k)} - s^{(k)}, k = 0, 1, 2, \dots,$

where $s^{(k)}$ is the solution of the system of linear equations

$J (f, x^{(k)}) \cdot s = f (x^{(k)}) .$

Jacobians are computed using the package ForwardDiff.jl. Plotting is done using the package Plots.jl.

18.9 μs

xxxxxxxxxx
 
begin
    using ForwardDiff
    using Plots
    plotly()
    using LinearAlgebra
end

18.6 s

Newton (generic function with 2 methods)

xxxxxxxxxx
 
function Newton(f::Function,J::Function,x::Vector{T},ϵ::Float64=1e-10) where T
    iter=0
    s=ones(T,size(x))
    ξ=x
    while norm(s)>ϵ && iter<100
        s=J(x)\f(x)
        ξ=x-s
        iter+=1
        x=ξ
    end
    ξ,iter
end

55.2 μs

Example 1

(Dennis and Schnabel, 1996) The solutions of the system

$\begin{aligned} 2 (x + y)^{2} + (x - y)^{2} - 8 & = 0 \\ 5 x^{2} + (y - 3)^{2} - 9 & = 0 \end{aligned}$

are the points $T_{1} = (1, 1)$ and $T_{2} \approx (- 1.18, 1.59)$ .

12.6 μs

f₁ (generic function with 1 method)

xxxxxxxxxx
 
# Vector function
f₁(x)=[2(x[1]+x[2])^2+(x[1]-x[2])^2-8,5*x[1]^2+(x[2]-3)^2-9]

29.9 μs

Float641

11.0

-3.0

 
f₁((1.0,2))

7.6 μs

Let us plot the functions and the contours in order to approximately locate the zeros:

9.5 μs

x
 
begin
    # Number of points
    m=101
    X=range(-2,stop=3,length=m)
    Y=range(-2,stop=3,length=m)
    # First applicate
    surface(X,Y,(x,y)->f₁([x,y])[1],xlabel="x",ylabel="y",colorbar=false)
    # Second applicate
    surface!(X,Y,(x,y)->f₁([x,y])[2],seriescolor=:blues)
end

232 ms

xxxxxxxxxx
 
begin
    # Locate the solutions using contour plot
    contour(X,Y,(x,y)->f₁([x,y])[1],contour_labels=true)
    contour!(X,Y,(x,y)->f₁([x,y])[2],contour_labels=true)
end

225 ms

xxxxxxxxxx
 
# Clearer picture
contour!(clims=(0,0.01),xlabel="x",ylabel="y",colorbar=:none)

459 μs

We see that the approximate zeros are $T_{1} = (- 1.2, 1.5)$ i $T_{2} = (1, 1)$ . Moreover, $T_{2}$ is exactly equal $(1, 1)$ (one iteration in the third computation). Furthermore, the method does not always converge (fourth computation).

22.7 μs

J₁ (generic function with 1 method)

xxxxxxxxxx
 
J₁(x)=ForwardDiff.jacobian(f₁,x)

19.3 μs

2×2 Array{Float64,2}:
 10.0  14.0
 10.0  -2.0

xxxxxxxxxx
 
# For example
J₁([1.0,2])

41.9 μs

11Float641

-1.18347

1.58684

21Float641

1.0

1.0

31Float641

1.0

1.0

41Float641

NaN

NaN

xxxxxxxxxx
 
Newton(f₁,J₁,[-1.0,0.0]), Newton(f₁,J₁,[0.5,1.1]), 
Newton(f₁,J₁,[1.0,1.0]), Newton(f₁,J₁,[0.0,0.0])

33.8 ms

Example 2

(Dennis and Schnabel, 1996) The solutions of the system

$\begin{aligned} x_{1}^{2} - x_{2}^{2} - 2 & = 0 \\ e^{x_{1} - 1} + x_{2}^{3} - 2 & = 0 \end{aligned}$

are the points $T_{1} = (1, 1)$ and $T_{2} \approx (- 0.71, 1.22)$ .

8.3 μs

xxxxxxxxxx
 
begin
    f₂(x)=[x[1]^2+x[2]^2-2,exp(x[1]-1)+x[2]^3-2]
    contour(X,Y,(x,y)->f₁([x,y])[1],contour_labels=true)
    contour!(X,Y,(x,y)->f₁([x,y])[2],contour_labels=true)
    contour!(clims=(0,0.01),xlabel="x",ylabel="y",colorbar=:none)
end

214 ms

11Float641

-0.713747

1.22089

21Float641

1.0

1.0

xxxxxxxxxx
 
begin
    J₂(x)=ForwardDiff.jacobian(f₂,x)
    Newton(f₂,J₂,[-1.0,1]), Newton(f₂,J₂,[0.8,1.2])
end

582 ms

Example 3

(Dennis and Schnabel, 1996) Solve $f (x) = 0$ , where

$f (x) = [\begin{matrix} x_{1} \\ x_{2}^{2} - x_{2} \\ e^{x_{3}} - 1 \end{matrix}] .$

The exact solutions are $T_{1} = (0, 0, 0)$ and $T_{2} = (0, - 1, 0)$ . We shall compute the solutions using several starting points.

28.6 μs

J₃ (generic function with 1 method)

 
begin
    f₃(x)=[x[1],x[2]^2+x[2],exp(x[3])-1]
    J₃(x)=ForwardDiff.jacobian(f₃,x)
end

31.1 μs

11Float641

0.0

0.0

0.0

21Float641

0.0

0.0

7.78375e-17

31Float641

0.0

0.0666667

NaN

41Float641

0.0

-1.0

0.0

xxxxxxxxxx
 
Newton(f₃,J₃,[-1.0,1.0,0.0]),Newton(f₃,J₃,[1.0,1,1]),
Newton(f₃,J₃,[-1.0,1,-10]),Newton(f₃,J₃,[0.5,-1.5,0])

717 ms

Example 4

(Rosenbrock parabolic valley) Given is the function

$f (x) = 100 (x_{2} - x_{1})^{2} + (1 - x_{1})^{2} .$

We want to find possible extremal points, thet is, we want to solve the equation

$g r a d f (x) = 0.$

15.5 μs

f₄ (generic function with 1 method)

 
f₄(x)=100(x[2]-x[1]^2)^2+(1-x[1])^2

27.3 μs

xxxxxxxxxx
 
# Plot the function using X and Y from Example 1
surface(X,Y,(x,y)->f₄([x,y]), seriescolor=:blues, xlabel="x", ylabel="y")

125 ms

xxxxxxxxxx
 
begin
    # This function is very demanding w.r.t. finding extremal points
    g₄(x)=ForwardDiff.gradient(f₄,collect(x))
    contour(X,Y,(x,y)->g₄([x,y])[1],contour_labels=true)
    contour!(X,Y,(x,y)->g₄([x,y])[2],contour_labels=true)
    contour!(clims=(-0.5,0.5),xlabel="x",ylabel="y",colorbar=:none)
end

272 ms

By observing contours, we conclude that this example is numerically demanding, while it is easy to see analytically that the only zero is $T = (1, 1)$ .

Here the vector function is given as the gradient of the scalar function $f$ , so the Jacobi matrix is computed using function FowardDiff.hessian() which approximates the matrix of second partial derivatives of the scalar function $f$ .

8.4 μs

1Float641

1.0

1.0

xxxxxxxxxx
 
Newton(g₄,x->ForwardDiff.hessian(f₄,x),[-1.0,2.0])

24.9 ms

Example 5

Let

$f (x) = \sum_{i = 1}^{11} (x_{3} \cdot \exp (- \frac{(t_{i} - x_{1})^{2}}{x_{2}}) - y_{i})^{2},$ where the pairs $(t_{i}, y_{i})$ are given by the table: $\begin{array}{cccccccccccc} i & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 \\ t_{i} & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 \\ y_{i} & 0.001 & .01 & .04 & .12 & .21 & .25 & .21 & .12 & .04 & .01 & .001 \end{array}$ We want to solve the equation $g r a d f (x) = 0.$ Unlike previous examples, where the condition number is $κ (J) = O (10)$

in Examples 1, 2 and 3 and

$κ (J) = O (1000)$

in Example 4, here

$κ (J) > O (10^{6})$

so the method is inaccurate and does not converge towrds the exact solution $x = (4.93, 2.62, 0.28)$ .

12.5 μs

f₅ (generic function with 1 method)

xxxxxxxxxx
 
begin
    t=collect(0:10)
    y=[0.001,0.01,0.04,0.12,0.21,0.25,0.21,0.12,0.04,0.01,0.001]
    f₅(x)=sum([( x[3]*exp(-((t[i]-x[1])^2/x[2]))-y[1])^2 for i=1:11])
end

3.3 ms

0.157753

2Float641

2.71553e-6

0.029986

1.13247

1.73704e5

xxxxxxxxxx
 
begin
    # Starting point is very near the solution
    x₀=[4.9,2.63,0.28]
    f₅(x₀)
    g₅(x)=ForwardDiff.gradient(f₅,x)
    J₅(x)=ForwardDiff.hessian(f₅,x)
    f₅(x₀), g₅(x₀), cond(J₅(x₀))
end

4.5 ms

1Float641

6.50244

0.00245085

-1.60888e-7

xxxxxxxxxx
 
x₅,iter₅=Newton(g₅,J₅,x₀,1e-8)

18.3 ms

Float641

1.52356e-51

2.04301e-49

-3.07371e-47

xxxxxxxxxx
 
g₅(x₅)

38.5 μs

1Float641

NaN

NaN

NaN

xxxxxxxxxx
 
Newton(g₅,J₅,[4.9,2.62,0.28],1e-8)

306 μs

Broyden's method

Given starting approximation $x_{0}$ and matrix $B_{0}$ , for each $k = 0, 1, 2, \dots$ , we compute:

$\begin{aligned} B_{k} \cdot s_{k} & = - f (x_{k}) (sustav) \\ x_{k + 1} & = x_{k} + s_{k} \\ y_{k} & = f (x_{k + 1}) - f (x_{k}) \\ B_{k + 1} & = B_{k} + \frac{(y_{k} - B_{k} s_{k}) s_{k}^{T}}{s_{k} \cdot s_{k}} \end{aligned}$

In this manner, we avoid the computation of the Jacobi matrix in each step. We can take $B_{0} = J (x_{0})$ , but also some other matrix.

30.7 μs

Broyden (generic function with 2 methods)

 
function Broyden(f::Function,B::Matrix,x::Vector{T},ϵ::Float64=1e-10) where T
    iter=0
    s=ones(T,length(x))
    ξ=x
    while norm(s)>ϵ && iter<100
        s=-(B\f(x))
        ξ=x+s
        y=f(ξ)-f(x)
        B=B+(y-B*s)*(s/(s⋅s))'
        x=ξ
        iter+=1
    end
    ξ,iter
end

58.3 μs

11Float641

-1.18347

1.58684

21Float641

1.0

1.0

xxxxxxxxxx
 
# Example 1
Broyden(f₁,J₁([-1.0,0.0]),[-1.0,0.0]), 
Broyden(f₁,J₁([1.0,1.5]),[1.0,1.5])

58.9 μs

11Float641

-0.0356391

-2.53166

21Float641

0.916022

-3.04547

31Float641

-1.18347

1.58684

xxxxxxxxxx
 
begin
    # Explain behaviour the method when we set B₀=I
    eye(n)=Matrix{Float64}(I,n,n)
    Broyden(f₁,eye(2),[-1.0,0.0]), Broyden(f₁,eye(2),[1.0,1.5]),
    Broyden(f₁,eye(2),[-1,1.5])
end

6.5 ms

11Float641

-0.713747

1.22089

21Float641

1.0

1.0

xxxxxxxxxx
 
begin
    # Example 2
    x0=[-1.0,1]
    x1=[0.8,1.2]
    Broyden(f₂,J₂([-1.0,1]),[-1.0,1]), Broyden(f₂,J₂([0.8,1.2]),[0.8,1.2])
end

70.8 μs

11Float641

0.0

5.96536e-26

0.0

21Float641

0.0

-1.0

0.0

xxxxxxxxxx
 
# Example 3
Broyden(f₃,J₃([-1.0,1,0]),[-1.0,1,0]), Broyden(f₃,J₃([0.5,-1.5,0]),[0.5,-1.5,0])

135 μs

11Float641

0.834064

0.695042

21Float641

1.0

1.0

31Float641

1.0

1.0

xxxxxxxxxx
 
# Example 4
Broyden(g₄,(x->ForwardDiff.hessian(f₄,x))([-1.0,2]),[-1.0,2]), # ali
Broyden(g₄,(x->ForwardDiff.hessian(f₄,x))([1,2.0]),[-1.0,2]),
Broyden(g₄,(x->ForwardDiff.hessian(f₄,x))([0.8,0.5]),[0.8,0.5])

13.0 ms

1Float641

18.9003

1.32399

0.0665517

xxxxxxxxxx
 
begin
    # Example 5
    x₆,iter₆=Broyden(g₅,(x->ForwardDiff.hessian(f₅,x))([4.9,2.6,0.2]), [4.9,2.6,0.2])
end

3.9 ms

Float641

1.85527e-29

-6.2359e-29

-2.07346e-29

xxxxxxxxxx
 
g₅(x₆)

45.9 μs

Davidon-Fletcher-Powell (DFP) method

DFP is an optimization method which finds extremal points of the function $F : R^{n} \to R$ , in which case $f (x) = g r a d F (x)$ .

Given initial approximation $x_{0}$ and matrix $H_{0}$ , for $k = 0, 1, 2, \dots$ , we compute:

$\begin{aligned} s_{k} & = - H_{k} f (x_{k}) \\ β_{k} & = \underset{β}{a r g m i n} F (x_{k} + β s_{k}) \\ s_{k} & = β_{k} s_{k} \\ x_{k + 1} & = x_{k} + s_{k} \\ y_{k} & = f (x_{k + 1}) - f (x_{k}) \\ H_{k + 1} & = H_{k} + \frac{s_{k} s_{k}^{T}}{y_{k} \cdot s_{k}} - \frac{H_{k} y_{k} y_{k}^{T} H_{k}}{y_{k} \cdot (H_{k} y_{k})} . \end{aligned}$

We can take $H_{0} = I$ . Notice that the iteration step does not require solving a system linear equations. Instead, all updates are performed using $O (n^{2})$ operations.

The one-dimensional minimzation along the line $x_{k} + β s_{k}$ is preformed by finding zeros of the directional derivative using bisection.

11.4 μs

Bisection (generic function with 2 methods)

xxxxxxxxxx
 
function Bisection(f::Function,a::T,b::T,ϵ::Float64=1e-10) where T
    fa=f(a)
    fb=f(b)
    x=T
    fx=T
    if fa*fb>zero(T)
        # return "Incorrect interval"
        if abs(fa)>abs(fb)
            return b,fb,0
        else
            return a,fa,0
        end
    end
    iter=0
    while b-a>ϵ && iter<1000
        x=(b+a)/2.0
        fx=f(x)
        if fa*fx<zero(T)
            b=x
            fb=fx
        else
            a=x
            fa=fx
        end
        iter+=1
        # @show x,fx
    end
    x,fx,iter
end

51.5 μs

DFP (generic function with 2 methods)

xxxxxxxxxx
 
function DFP(f::Function,H::Matrix,x::Vector{T},ϵ::Float64=1e-10) where T
    iter=0
    s=ones(T,length(x))
    ξ=x
    while norm(s)>ϵ && iter<50
        s=-H*f(x)
        s0=s/norm(s)
        F(ζ)=f(x+ζ*s)⋅s0
        β,fx,iterb=Bisection(F,0.0,1.0,10*eps())
        s*=β
        ξ=x+s
        y=f(ξ)-f(x)
        z=H*y
        H=H+(s/(y⋅s))*s'-(z/(y⋅z))*z'
        x=ξ
        iter+=1
    end
    ξ,iter
end

105 μs

Example 6

Let us find extremal points of the function

$f (x, y) = (x + 2 y - 7)^{2} + (2 x + y - 5)^{2} .$

The function has minimum at $(1, 3)$ .

13.2 μs

f₆ (generic function with 1 method)

xxxxxxxxxx
 
f₆(x) = (x[1] + 2*x[2]-7)^2 + (2*x[1] + x[2]-5)^2

49.5 μs

xxxxxxxxxx
 
f₆([1,2])

1.0 μs

g₆ (generic function with 1 method)

xxxxxxxxxx
 
g₆(x)=ForwardDiff.gradient(f₆,x)

21.2 μs

1Float641

1.0

3.0

xxxxxxxxxx
 
DFP(g₆,eye(2),[0.8,2.7],eps())

47.0 ms

1Float641

1.0

1.0

xxxxxxxxxx
 
# Example 4
DFP(g₄,eye(2),[0.9,1.1])

293 μs

1Float641

4.89611

46.1979

0.00158913

xxxxxxxxxx
 
# Example 5
DFP(g₅,eye(3),[4.9,2.6,0.2])

2.3 ms

Broyden-Fletcher-Goldfarb-Schano (BFGS) method

BFGS is an optimization method for finding extremal points of the function $F : R^{n} \to R$ , in which case $f (x) = g r a d F (x)$ .

The method is similar to the DFP method, with somewhat better convergence properties.

Let $F : R^{n} \to R$ , whose minmum we seek, and let $f (x) = g r a d F (x)$ .

Given initial approximation $x_{0}$ and matrix $H_{0}$ , for $k = 0, 1, 2, \dots$ , we compute:

$\begin{aligned} s_{k} & = - H_{k} f (x_{k}) \\ β_{k} & = a r g m i n F (x_{k} + β_{k} s_{k}) \\ s_{k} & = β_{k} s_{k} \\ x_{k + 1} & = x_{k} + s_{k} \\ y_{k} & = f (x_{k + 1}) - f (x_{k}) \\ H_{k + 1} & = (I - \frac{s_{k} y_{k}^{T}}{y_{k} \cdot s_{k}}) H_{k} (I - \frac{y_{k} s_{k}^{T}}{y_{k} \cdot s_{k}}) + \frac{s_{k} s_{k}^{T}}{y_{k} \cdot s_{k}} . \end{aligned}$

We can take $H_{0} = I$ . Notice that the iteration step does not require solving a system linear equations. Instead, all updates are performed using $O (n^{2})$ operations.

The one-dimensional minimzation along the line $x_{k} + β s_{k}$ is preformed by finding zeros of the directional derivative using bisection.

22.5 μs

BFGS (generic function with 2 methods)

xxxxxxxxxx
 
function BFGS(f::Function,H::Matrix,x::Vector{T},ϵ::Float64=1e-10) where T
    iter=0
    s=ones(T,length(x))
    ξ=x
    while norm(s)>ϵ && iter<50
        s=-H*f(x)
        s0=s/norm(s)
        F(ζ)=f(x+ζ*s)⋅s0
        β,fx,iterb=Bisection(F,0.0,1.0,10*eps())
        s*=β
        ξ=x+s
        y=f(ξ)-f(x)
        z=H*y
        α=y⋅s
        s1=s/α
        H=H-s1*z'-z*s1'+s1*(y⋅z)*s1'+s1*s'
        x=ξ
        iter+=1
    end
    ξ,iter
end

104 μs

1Float641

1.0

3.0

xxxxxxxxxx
 
BFGS(g₆,eye(2),[0.8,2.7],eps())

50.8 ms

1Float641

1.0

1.0

xxxxxxxxxx
 
# Example 4
BFGS(g₄,eye(2),[0.9,1.1])

264 μs

1Float641

2.08908

31226.6

0.00100056

xxxxxxxxxx
 
# Example 5
BFGS(g₅,eye(3),[4.9,2.6,0.2])

2.1 ms

Julia packages

Previous programs are simple illustrative implementations of the mentioned algorithms. Julia package NLsolve.jl is used for solving systems of linear equations and the package Optim.jl is used for nonlinear optimization.

13.1 μs

xxxxxxxxxx
 
using NLsolve

9.7 s

f₁! (generic function with 1 method)

xxxxxxxxxx
 
# Example 1
function f₁!(fvec,x)
    fvec[1] = 2(x[1]+x[2])^2+(x[1]-x[2])^2-8
    fvec[2] = 5*x[1]^2+(x[2]-3)^2-9
end

45.3 μs

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [-1.0, 0.0]
 * Zero: [-1.1834670032425283, 1.5868371427230779]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6

xxxxxxxxxx
 
nlsolve(f₁!,[-1.0,0])

178 ms

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.5, 1.1]
 * Zero: [1.0000000000000002, 1.0000000000000002]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6

xxxxxxxxxx
 
nlsolve(f₁!,[0.5,1.1])

40.7 μs

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [-1.0, 1.0]
 * Zero: [-0.7137474114758742, 1.2208868222037403]
 * Inf-norm of residuals: 0.000000
 * Iterations: 4
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 5
 * Jacobian Calls (df/dx): 5

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.8, 1.2]
 * Zero: [0.9999999999940328, 1.0000000000115203]
 * Inf-norm of residuals: 0.000000
 * Iterations: 4
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 5
 * Jacobian Calls (df/dx): 5

xxxxxxxxxx
 
begin
    # Example 2
    function f₂!(fvec,x)
        fvec[1] = x[1]^2+x[2]^2-2
        fvec[2] = exp(x[1]-1)+x[2]^3-2
    end
    nlsolve(f₂!,[-1.0,1]), nlsolve(f₂!,[0.8,1.2])
end

168 ms

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [-1.0, 1.0, 0.0]
 * Zero: [0.0, 2.3283064364709337e-10, 0.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.0, 1.0, 1.0]
 * Zero: [0.0, 2.3283064364709337e-10, 1.223235687436471e-12]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [-1.0, 1.0, -10.0]
 * Zero: [0.0, 4.9978109751571114e-11, 8.131707812509303e-17]
 * Inf-norm of residuals: 0.000000
 * Iterations: 20
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 21
 * Jacobian Calls (df/dx): 8

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.5, -1.5, 0.0]
 * Zero: [0.0, -1.0000000000000007, 0.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 6

xxxxxxxxxx
 
begin
    # Example 3
    function f₃!(fvec,x)
        fvec[1] = x[1]
        fvec[2] = x[2]^2+x[2]
        fvec[3] = exp(x[3])-1
    end
    nlsolve(f₃!,[-1.0,1.0,0.0]), nlsolve(f₃!,[1.0,1,1]),
    nlsolve(f₃!,[-1.0,1,-10]), nlsolve(f₃!,[0.5,-1.5,0])
end

153 ms

xxxxxxxxxx
 
using Optim

8.4 s

 * Status: success

 * Candidate solution
    Final objective value:     5.375030e-17

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 5.13e-09 ≰ 0.0e+00
    |x - x'|/|x'|          = 5.13e-09 ≰ 0.0e+00
    |f(x) - f(x')|         = 9.67e-17 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.80e+00 ≰ 0.0e+00
    |g(x)|                 = 2.10e-11 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    35
    f(x) calls:    102
    ∇f(x) calls:   102

xxxxxxxxxx
 
# Example 4
optimize(f₄,[-1.0,2],Optim.BFGS())

382 μs

 * Status: success

 * Candidate solution
    Final objective value:     3.572548e-12

 * Found with
    Algorithm:     BFGS

 * Convergence measures
    |x - x'|               = 1.87e+05 ≰ 0.0e+00
    |x - x'|/|x'|          = 8.70e-02 ≰ 0.0e+00
    |f(x) - f(x')|         = 6.03e-16 ≰ 0.0e+00
    |f(x) - f(x')|/|f(x')| = 1.69e-04 ≰ 0.0e+00
    |g(x)|                 = 9.69e-09 ≤ 1.0e-08

 * Work counters
    Seconds run:   0  (vs limit Inf)
    Iterations:    41
    f(x) calls:    137
    ∇f(x) calls:   137

xxxxxxxxxx
 
# Example 5 - again there is no convergence
optimize(f₅,[4.9,2.6,0.2],Optim.BFGS())

4.3 ms

Enter cell code...xxxxxxxxxx
 
​

200 ns